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Abstract 

Here we present/implement an algorithm to find Liouvillian first integrals of 
dynamical systems in the plane. In [1], we have introduced the basis for the present 
implementation. The particular form of such systems allows reducing it to a single 
rational first order ordinary differential equation (rational first order ODE). We 
present a set of software routines in Maple 10 for solving rational first order ODEs. 
The package present commands permitting research incursions of some algebraic 
properties of the system that is being studied. 
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PROGRAM SUMMARY 



Title of the software package: Lsolver. 
Catalogue number: (supplied by Elsevier) 

Software obtainable from: CPC Program Library, Queen's University of Belfast, N. Ireland. 
Licensing provisions: none 

Operating systems under which the program has been tested: Windows ME, Windows XP. 
Programming languages used: Maple 10 

Memory required to execute with typical data: 128 Megabytes. 

No. of lines in distributed program, including On-Line Help, etc.: 900. 

Keywords: Liouvillian functions, first integrals, dynamical systems in the plane, first order 
ordinary differential equations, computer algebra, Prelle-Singer (PS). 

Nature of mathematical problem 

Solution of rational first order ordinary differential equations. 
Methods of solution 

The method of solution is based on a Darboux/PS type approach. 
Restrictions concerning the complexity of the problem 

If the integrating factor for the rational first order ODE under consideration presents Dar- 
boux Polynomials of high degree ( > 3 ) in the dependent and independent variables, the 
package may spend an unpractical amount of time to obtain the solution. 

Typical running time 

This depends strongly on the ODE, but usually under 2 seconds. 
Unusual features of the program 

Our implementation not only search for Liouvillian conserved quantities, but can also be 
used as a research tool that allows the user to follow all the steps of the Darboux procedure 
(for example, the algebraic invariants curves and associated cofactors can be calculated). In 
addition, since our package is based in recent theoretical developments [1], it can successfully 
solve a class of rational first order ODEs that were not solved by some of the best-known 
ODE solvers available. 



2 



LONG WRITE-UP 



1 Introduction 

The problem of searching for conserved quantities (first integrals) on physical sys- 
tems is not new. In the last thirty years there has been a great increasing of interest 
in methods that follow the lines drawn by Darboux [2] in the second half of 19 th 
century combined with topics that comes from differential algebra (developed by 
Ritt [3] and Kolchin[4] in the middle of the 20 th century). 
Particularly, the polynomial systems on the plane, 

Tt = N {x ,y) 

f = M(x,y) (1) 

where M and N are polynomials in (x,y), have received a great deal of attention. 
We can then notice that there is an equivalence in finding first integrals of this 
system with the solving of a rational first order differential equation of the form: 

, = dy = M(x,y) 
V dx N(x,y)' [ ' 

If I(x, y) is a first integral of the system (1) then I(x, y) = K, where K is a constant, 
is a general solution of the first order differential equation (2). 

A big step forward in constructing an algorithm for solving first order ODEs an- 
alytically was taken in a seminal paper by Prelle and Singer (PS) [5] on autonomous 
systems of ODEs. Prelle and Singer's problem is equivalent to asking when a ra- 
tional first order ODE of the form (2) has an elementary solution (a solution which 
can be written in terms of a combination of polynomials, logarithms, exponentials 
and radicals). Prelle and Singer were not exactly able to construct an algorithm 
for solving their problem, since they were not able to define a degree bound for the 
polynomials which might enter into the solution. Though important from a theo- 
retical point of view, any computer implementation of the PS method will have a 
practical degree bound imposed by the processing speed and/or memory needed to 
handle the ever-more complex equations. With this in mind it is possible to say that 
Prelle and Singer's original method is almost an algorithm, awaiting a theoretical 
degree bound to turn it algorithmic. The approach of Prelle and Singer [5] was so 
interesting that it has motivated many extensions to its main ideas [6, 7, 8, 9, 10, 11] 

The purpose of this paper is to describe an implementation, in Maple, of a 
Darboux- type procedure (presented on [1]) that extends the applicability of the 
PS-approach to deal with rational first order ODEs (2) with Liouvillian 2 solutions. 

The paper is organized as follows: in section 2, we present a short theoretical 
introduction to the PS approach; in the following section, we introduce the main 
ideas of our extension [1]; in section 4, we present examples of the application of 
our method; in section 5, we present the Lsolver package; section 6 talk about some 
interesting features and performance of our package. 



2 For a formal definition of Elementary and Liouvillian functions, please see [12]. 
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2 The State of the art 



Despite its usefulness in solving rational first order ODEs, the Prelle-Singer proce- 
dure is not very well known outside mathematical circles, and so we present a brief 
overview of the main ideas of the original PS-procedure and some of its extensions. 



2.1 The Prelle-Singer Method 

Consider the class of rational first order ODEs which can be written as 

, = dy = M(x,y) 
V dx N(x,y) { ' 

where M(x, y) and N(x, y) are polynomials with coefficients in the complex field C. 

In [5], Prelle and Singer proved that, if an elementary first integral of (3) exists, 
it is possible to find an algebraic integrating factor R such that R n is a rational 
function of (x, y) for some integer n. If one knows R, then the ODE can be solved 
by quadrature. 

We have that: 

a(|A0 + a(|M) =0 

ox dy 

From (4) we see that 

dR dN dR dM 

N 7T + R ir + M ir + R ^ = - 5 

ox Ox oy Oy 

Thus 

D[R] __ /dN dM" 
R \ dx dy 

where 

and 

i 

where fi are irreducible polynomials and n-i are non-zero rational numbers. From 
(7) and (8), we have 



+ ST , (6) 



d[r] D[Uifn ^ifr^iDmn^ii 
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R Uifk k U k fk k 

i fi' i fi 

From (6), plus the fact that M and N are polynomials, we conclude that D[R]/R 
is a polynomial. Therefore, from (9), we see that fi\D[fi] (i.e., fi is a divisor of 
D[fi]). 

We now have a criterion for choosing the possible fi (build all the possible 
divisors of D[fi]) and, by using (6) and (9), we have 

If we manage to solve (10) and thereby find nj, we know the integrating factor 
for the ODE and the problem is reduced to a quadrature. 
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2.2 Extensions to the Prelle and Singer method 



The original PS-procedure deals with rational first order ODEs with elementary 
solutions. Extensions were produced, both on the theoretical basis and on producing 
new algorithms, in order to allow tackling more general ODEs. 

Using the results presented on [13], it is possible to shown [9, 7, 8] that, if 
the solution for a rational first order ODE of form (3) can be written in terms of 
Liouvillian functions, then it presents an integrating factor is of the form: 

n 

R = e r (x,y) -Q y yi = e P(x,y)/Q(x,y) y ^ (n) 

i=l 

where ro = P/Q (P and Q are polynomials in (x, y)) is a rational function of (x, y), 
D[ro] is a polynomial, S = Ylj v//\ the V{ are Darboux polynomials (eigenpolynomi- 
als) of the D operator and q are constants. 

Note that this result is the analogous of the previously mentioned result due 
to Prelle-Singer but for a more general case, namely, now we are dealing with 
Liouvillian solutions (please note that the elementary solutions, dealt with by Prelle 
and Singer, are also Liouvillian but are a restriction to the general Liouvillian case). 

In [7, 8], we introduced a method that, using the knowledge embodied by (11), 
solves the ODE for the cases where ro is either f(x), g{y) or f(x) + g(y), where / 
and g are rational functions. 

In the next section, we will go a further step by introducing a new method that 
allows for the general case for the rational function tq. 



3 Introducing the new Method 

In what follows, we will refine the result (11) by uncovering information about the 
structure of Q (ro = P/Q)- This further knowledge will then be converted into 
something that will allow us to build a solving method. 
Using (11) into D[R]/R = —(d x N + d v M), we get: 

+ ^ = -(d y M + d x N), (12) 
leading to 

om^m +Eej £M = - ia ^ +vn . (13 , 

Remembering that the Vj are Darboux polynomials of the D operator, we can 
write D[vj] = Xj Vj, where the Xj are polynomials in (x,y) (called cofactors) associ- 
ated with the vj. Thus, we can write (13), multiplying both sides by Q 2 , as: 

Q D[P] - P D[Q] + Q 2 ]T Cj Xj = -Q 2 (d y M + d x N) . (14) 

j 

The above equation would be the analogous (in our method) to equation (6) (in 
the PS method). We call the attention to the fact that we are trying to deal with a 
much wider class of ODEs. This is embodied by the presence of the term e p ^ on 
the integrating factor. So far, we only know that P and Q are polynomials thus, it 
is easy to see that, attempting to solve (14), actually means solving a third degree 
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equation for the unknowns: Cj and the coefficients defining the polynomials P and 
Q. So, our case here seems much poorer than the PS original one as a practical 
approach to solving rational first order ODEs (albeit dealing with more general 
ones) . 

Can we do something to improve that situation? Actually, we will show that we 
can by presenting a refinement to our knowledge about the general structure for the 
integrating factor for ODEs of the form (3). 



Theorem: Let the exponent tq be expressed as P(x,y)/Q(x,y), where P and Q 
are polynomials in (x,y), with no common factor. Then we have that Q\D[Q] (i.e., 
D[Q]/Q is a polynomial in (x,y)). 

Proof: Since D[ro] is polynomial (see [8]), we can write it as: 



D[r ] = D 



QD[P]-PD[Q] =U (15) 



Q 

where n is polynomial in (x,y). Multiplying (15) by Q, one obtains: 

D[P] _^m =IlQ (16) 

Since D is a linear differential operator, with polynomial coefficients, and P is 
polynomial, D[P] is also polynomial. Therefore, since HQ is also polynomial, we 
may conclude that PD q^ is polynomial either. Since, by hypothesis, P and Q have 
no common factor, we can infer that Q\D[Q] (i.e., D[Q]/Q is a polynomial), as we 
wanted to demonstrate 3 . 



This result, in turn, leads to the following corollary: 

Corollary: We can write Q as Yli=iQi( x ^y) mi > where the qi are irreducible 
independent Darboux polynomials of the D operator and the rrii are positive integers. 



Proof: Q is a polynomial so it can be written as T\i=i qi(x , y) mi , where the qi 
are independent irreducible polynomials and the m^'s are positive integers. So, we 
can write: 

^ = (17) 

H i Qi 

from the theorem just proved, we know that the left-hand side of (17) is polynomial. 
So the right-hand side is also polynomial implying that <Zi|£%i] (i.e., D[qi]/qi is a 
polynomial), as we wanted to show. 



Using these results, we can write Q = JJi {vqij m \ where the v q j are independent 
irreducible Darboux polynomials of the D operator, and so: 

as =I>( ss! =I>v (18) 

Q i v qi 



3 An analogous result was demonstrated (by other means) on [10] 
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where the are the cofactors associated with the v^. Using this into (14) and 
re-arranging we get: 

D[P] - P J2 mi A gi = - J] (vqi) nH I J2 c 3 ^ + d y M + 8 x n\ . (19) 

i i \ j ) 

This equation will prove to be very important to our method. So, let us now do 
some analyzing of its structure: What is the advantage of (19) in comparison with 
(14)? 

In words, now we know the structure of Q, i.e., the building blocks of Q. In (14), 
we had to set the degree for Vi (analogously, in the PS method, we had to set the 
degree for the fi) and to set the degree for the P and Q polynomials (there is no 
analogous in the PS method). After that, as already pointed out, we are left to solve 
a third degree equation on the unknowns. Now, with (19), we have to do the same 
setting of degrees but, due to our improved knowledge about Q, after setting these 
degrees, we are left with a finite set of possibilities for the integers {rrii}. Thus 
converting solving (19) to solving a system of linear equations for the unknowns Cj 
and the coefficients defining P (one system of linear equations for each possible set 
of integers {rrii}). 

So, after setting the degrees we mentioned, in essence, equation (19) is as man- 
ageable as equation (6) of the PS method. This is the basis of a new method to 
deal with rational first order ODEs that can be summarized as follows: 

In equation (19) (the cornerstone of our method), we see terms involving the 
Darboux polynomials (and/or the associated cofactors) of the D operator. So, the 
first job at hand is to determine those. In order to do that we have to choose a 
degree (as in the PS method). For simplicity sake, let us start with the most simple 
possibility, namely, find the Darboux polynomials of degree 1 of the D operator 
and the associated cofactors. To be able to solve equation (19), we have to set the 
degree for the Q and P polynomials. Remembering that we know the structure 
for Q, setting the degree for Q allows us (see above) to find all possibilities for 
{rrii}. After this, we then try to solve (19). If we solve it, we would have found the 
integrating factor we were looking for. If we can not solve it, we have to change our 
settings and try again. The way we will change the settings 4 is a matter of choice 
and, in principle, it is difficult to see which one will be the best since this depends 
strongly on the ODE itself. The point is that, no matter which is that choice, the 
method is contained in the sense that, even if we have to use high values for the 
degree for Q and P, our task will be to solve linear equations (sure, it can be a 
great number of equations, but always linear ones). 

Of course, this procedure can go on indefinitely, but we will be covering all the 
possibilities and we may hope that we will find a solution within our lifetime. On 
a brighter tone, all the examples we have come across are solved with degree for 
the Darboux polynomials of the D operator, Q and P that allow for a practical 
realization of the method. 

4 For example, should we increase the degree of P without increasing the one for Q? Up to which 
value? When should we increase the degree for the Darboux polynomials of the D operator? etc. 
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4 The Inner Workings of the Method 



In this section, we are going to present examples of application of the method. We 
will present the calculations in a long version in order for the reader to have a deeper 
understanding of the inner workings of the method. 

First, in order to illustrate the method just presented, we are going to start with 
a simple ODE. This example was chosen due to its simplicity (so it is a good stating 
point for the introducing of the method) and for the fact that, even being simple, 
this ODE is not solved by the Maple powerful solver dsolve. 

Example 1: 

Consider the following ODE: 

dy = (x+l)y 
dx x — xy — y 2 + x 2 

For this equation, we have M = (x + 1) y and N = x — xy — y 2 + x 2 leading to: 

D = N8 X + Md y = (x-xy-y 2 + x 2 )d x + (x + 1) yd y (21) 

Up to degree 1, we have that the Darboux polynomials (with the associated cofac- 
tors) for this operator are: 

• v-i = y, Ai = x+ 1, 

• V2 = x + y, A 2 = 1 + x - y. 

Then we have to choose the degree for the polynomial Q. Starting with degree 1, 
since v\ and i> 2 are of degree 1, the only possible values for rrii are {mi = 1, to 2 = 0} 
and {mi = 0, to 2 = 1}. Starting with degree 1 for P, we have P = a\ + a 2 x + 03 y 
and equation (19) leads to: 



y ml {x + y) m ^{nl {x + 1) + n2 (1 + x - y) + 3 x + 2 - y) + 
(x - xy - y 2 + x 2 ^ a2 + (x + 1) ya3 - 
(al +a2x + a3y) {ml (x + 1) + m2 (1 + x - y)) = 0. (22) 

As we shall see, with the values m\ = 1, m 2 = it is possible to find a solution. 
Substituting {mi = l,m 2 = 0} into (22) we get: 

y (nl (x + 1) + n2 (1 + x — y) + 3 x + 2 — y) + (^x — xy — y 2 + x 2 ^j a2 + 

(x + 1) ya3 -{al + a2x + a3 y) (x + l) = 0. (23) 

In order to solve the above equation, the coefficients for different powers of {x, y) 
have to be zero. Thus leading to the following system of linear equations: 



nl + n2 + 2 = 
-n2 - a2 - 1 = 
-al = 
nl + n2 + 3 - a2 = 



(24) 
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Leading to the solution for the coefficients: 

al = 0, a3 = a3, nl = 0, ra2 = -2, a2 = 1 (25) 
So, the integrating factor for this ODE becomes (choosing 03 = 0): 

P x/v 

R = — — , (26) 
(x + y) 

Then the solution is: 

C = v(-} + v)f ,V _ e -i K(1 , _£±») (27 ) 
(x + y) ' y 

We are going to present now an example extracted from the book by Kamke 
[14]. This choice for the second example was made aiming to show that, even for 
involved cases, our method is contained and manageable. 

Example 2: 

Consider the following ODE (1.169 from the book by Kamke): 

(ax + bf ^ + (ax + b)y 3 + cy 2 (28) 
ax 

For this equation, we have M = — ((ax + b)y 3 + cy 2 ) and N = (ax + 6) 2 leading 

to: 

D = Nd x + Md y = (ax + b) 2 8 X - ((ax + b)y 3 + cy 2 ) d y (29) 

For this equation, up to degree 1, we have that the Darboux polynomials (with 
the associated cofactors) for this operator are: 

• vi = y, Ai = -yc- by 2 - axy 2 , 

• V2 = (ax + b)/a, \2 = ab + a 2 x. 

Then we have to choose the degree for the polynomial Q. For this particular 
example, we will see that we need to use degree 4 for both Q and P. Since v\ 
and V2 are of degree 1, the possible values for rrii are {m\ = 4, rri2 = 0}, {m\ = 
3, m 2 = 1}, {mi = 2,m 2 = 2}, {mi = 1,777.2 = 3} and {ttt-i = 0,tt7 2 = 4}. Letting 
P = 01 y A + a 2 x 2 y + a 3 t/ 2 + a 4 x y 2 + a 5 x y + a 6 x 4 + a 7 x 3 + a 8 x 2 + a 9 y 3 + ai x 3 y + 
a n xy 3 + ai2 x + ai3x 2 y 2 + 014 y + 015, equation (19) becomes: 

y mi ^ ax + ^_ yc _ by 2 _ axy 2^ + n2 (ab + a 2 x^j 

-2 y (axy + by + c) -y 2 (ax + 6) + 2 a 2 x + 2 aft) 
+ (a 2 x 2 + 2axb + 6 2 ) (2 xy + a4 y 2 + a5 y+ 
4a6x 3 + 3a7x 2 + 2a8x + 3al0x 2 y + ally 3 + al2 + 2 al3 xy 2 ) 
-y 2 (axy + by + c) (a al y 3 + a2 x 2 + 2 a5 y + 2 xy + a5 x + 3 a9 y 2 



+ al0 x 3 + 3 all xy 2 + 2al3 x 2 y + ai^) - (al y 4 + x 2 y + a3 y 2 + a4 xy 2 
+ a5 xy + a6 x 4 + a7x 3 + a8 x 2 + a9 y 3 + alO x 3 y + aii xy 3 + aiS x 
+ a!5 x 2 y 2 + al4y+ al5^j (ml (-yc - by 2 - axy 2 ) + m2 (ab + a 2 x)) = (30) 
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As we shall see, with the values mi = 2, m,2 = 2, it is possible to find a solution. 
Substituting those values for m\ and m2 into (30) we get: 

[Ab 2 a6 + 6aba7 + 2a 2 a8 = 0,3b 2 a7 + a 2 al2 + Aaba8 = 0, 

2aba5 + al2c+n2a 2 + 2b 2 a2 + 2a 2 = 0,-3cal - 2ba9 = 0, 
b 2 al2 = 0,-nl c+ al5b + b 2 a4 - 2c = 0,2a 2 a2 + a7c + 6abal0 = 0, 
8aba6 + 3a 2 a7 = 0,2b 2 a8 + 2abal2 = 0,Aa 2 a6 = 0, 
-3 6 -nl b + b 2 all - ca3 = 0,a 2 a5 + 3b 2 al0 + a8c + 4aba2 = 0, 
bal2 + al5 a + 2b 2 al3 + 2aba4 = 0, -nl a-3a-ca4 +2aball = 0, 
ba8 + a 2 a4 + aal2 + Aabal3 = 0,n2 ab + 2ab + b 2 a5 + al5 c = 0, 
a 2 all - calS = 0,-ba3 -2ca9 = 0,a6 a = 0,3a 2 alO + a6c = 0, 
-2ball - 2aa9 = 0,-3aal =0,-3bal =0,-aal3 = 0,-2aall =0, 

-bal3 - aa4 = 0, -2call - aa3 - ba4 = 0, 

aa8 +ba7 + 2a 2 al3 = 0,aa7 + ba6 = o} (31) 

Leading to the solution for the coefficients: 

{a6 = 0, all = 0, alO = 0, al = 0, n2 = -1, a2 = 0, a7 = 0, a9 = 0, 

i , c 2 -2ab 2 al3 n cb 

a3 = -1/2 = , nl = -3, al4 = 

a 6 a A 

b 2 c 
al5 = -1/2 — , a5 = — , al2 = -b, a4 = 2 
a a 



bal3 



} al3 = al3,a8 = -a/2 J (32) 



So, the integrating factor for this ODE becomes (choosing 013 = 0): 

(^yc-\-a^ x-\-ab^j 
g 2 ay 2 (ax + b) 2 

R = — ^ — (33) 



y 3 (ax + b) 



Then the solution is: 



c = ( ^y+^y+^li^ *_i dx (34) 

J y(ax + b) 

From the example above, we can see that, even in cases where Q and P have 
high degrees, our method convert the problem into a simple linear algebraic system. 
More specifically, in this example, the equations (31). Although these equations 
look frightening, they are solvable even by hand (if one is brave enough). In this 
computational implementation, the problem is trivially solved. 

Let us then introduce our package implementing the above explained ideas and 
making the algorithm available for easy usage. 
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5 The Lsolver package 
Summary 

A brief review of the commands of the package is as follows: 5 

• Lsolve solves ODEs, using our Darboux-type approach. It deals only with 
rational first order ODEs. 

• IntFact returns an integrating factor for the rational first order ODE. 

• LDop constructs the D operator associated to the rational first order ODE. 

• Darboux returns the DPs (Darboux polynomials) and cofactors associated with 
the D operator. 

Description 

A complete description of the Lsolver package's commands can also be found in the 
on-line help. 

5.1 Command name: Lsolve 

Feature: This command applies our Darboux-type procedure to solve ODEs. 

Calling sequence 6 : 

> Lsolve (ODE, extra_args) ; 

Parameters: 

ODE - a first order ordinary differential equation. 

extra_args - Deg = <bound> - a positive integer equal to the maximum degree of 
the darboux polynomials. The default value is 1. 

- Deg_Q = <bound_Q> - a positive integer equal to the maximum degree of 
the Q polynomial; The default value is 5. 

- Deg_P = <bound_P> - a positive integer equal to the maximum degree 
of the P polynomial; The default value is 5. 

Synopsis: 

The Lsolve command is a part of the Lsolver package, which is designed to solve 
rational first order ODEs. 

The arguments 

The first argument of Lsolve is the ODE. The extra arguments are all positive 
integers that specify the polynomial degrees we are going to use in the search for 

5 This subsection and the next one may contain some information already presented in the previous 
sections; this is necessary to produce a self-contained description of the package. 
6 In what follows, the input can be recognized by the Maple prompt >. 
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the integrating factor. We would like to remind the reader that our Darboux-type 
procedure is a semi-decision one, in other words, if the integrating factor has not be 
found for these degrees, it does not exists at this level. But, it could as well exist 
for higher degree, we have to keep looking. 

5.2 Command name: IntFact 

Feature: Calculates an integrating factor for the ODE. 
Calling sequence: 

> IntFact (ODE, extra_args) ; 

Parameters: 

The parameters ODE and extra_args have the same meaning as explained above . 
The command IntFact admits all the extra_args as the command Lsolve. 

Synopsis: 

Actually this command embodies the heart of the package since, after finding 
the integrating factor, the solution is found via a quadrature. 

5.3 Command name: Darboux 

Feature: Calculates the Darboux polynomials, up to a certain degree, of the D 
operator related to the ODE in question. The default is to use degree=l for the 
Darboux polynomials to be sought. 

Calling sequence: 

> Darboux ( ODE, extra_args) ; 

Parameters: 

The command Darboux admits as extra_args only Deg = <bound>, where <bound> 
is a positive integer equal to the maximum degree of the darboux polynomials. 

Synopsis: 

This command is very important to the method since it determines the building 
blocks for the integrating factor. 

5.4 Command name: LDop 

Feature: Returns the D operator (for the ODE) that is an essential ingredient in 
applying our procedure. 

Calling sequence: 

> Dop(ODE) ; 

Synopsis: 

This command returns the D operator associated with the ODE. 
7 In what follows, these parameters will always have this meaning. 



12 



5.5 Example of the usage of the package commands 

Let us use some rational first order ODEs to show the commands in action. 
Consider the following ODE: 

dy y{l+x) 



(35) 

ax x + ar — y z 

Usually, people only want to find the solution for the ODE. So, that would 
simply require the input lines: 

> eq := diff(y(x),x) = y(x)*(l+x)/(x+x~2-y(x)~2) ; 

> Lsolve (eq) ; 

The solution will be 

-zvW2er/ ( l l 2l ^ x \ + 2 e 1/2 l£^ y ( x ) = _C1 (36) 

V y( x > J 

Sometimes the user is interested on having a look at just the integrating factor, 
so he/she could type: 

> IntFact (eq) ; 

For the ODE in the form 



dy y(l + x) 



dx x + x 2 — y 2 
the integrating factor will be 



exp 



X 2 

■ill 1 



y 2 



(37) 



(38) 



Please notice that the program informs the user the format in which it is considering 
the ODE to be in, that is important since it may be the case where the input ODE 
gets transformed by the simplification routines and is regarded in a different shape 
(of course it is still the same ODE) and that will affect the integrating factor. 

If one wants to find the Darboux polynomials (up to a certain degree), one can 
use the following input: 

> Darboux (eq) ; 

Below, you can find a list containing two lists: 
[[Darboux polynomials] , [co-factors]] 



[&/],[!+*]] (39) 

Please remember that the default value for the degree of the Darboux polynomials 
is 1. 

Finally, for this example, let us calculate the D operator: 
> LDop(eq) : 
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w -► {x + x 2 -y 2 )(^ w) + y{l + x){-jj- w) (40) 
Let us deal with a more involved case. Consider the following ODE. 

dy = -l + x + y + 3y 2 
dx 2(2 x + y + xy + y 2 — y 3 ) 
Again, using the input lines just introduced above: 

> eq := dif f (y (x) ,x) = l/2*(-l+x+y(x)+3*y(x) ~2)/(2*x+y(x)+x*y(x)+ 
y(x)~2-y(x)~3) ; 

> Lsolve(eq) ; 

We could not find an integrating factor 

What happened? The reason for the failure in solving the equation above could 
be due to the fact that there is no Darboux polynomials with degree equals to 1 
and, since the default for the command assumes that value, we have to increase this 
value to check that possibility. This is done via the following input line: 

> Lsolve(eq,Deg=2) ; 

The solution will be 



/ i I A -l+4x+4y(x) 1 _1 I O / \ \ 

2Jx + (y(x)) 2 e 1 *Hv(*))* erf 1/2 + y ^ = _C1 (42) 

V vW(y(*))V 

Let us have a look at the other commands: 

> IntFact (eq,Deg=2) ; 

For the ODE in the form 

dy -l + x + y + 3y 2 



(43) 



dx Ax + 2y + 2yx + 2y 2 -2y 3 
the integrating factor will be 

-1/4+x + y 9 \-3/2 

e x +v 2 [x + y z \ (44) 
The other two commands, for this case, produce: 

> Darboux (eq,Deg=2) ; 

Below, you can find a list containing two lists: 
[[Darboux polynomials] , [co-factors]] 

[[x + y 2 ],[4 + 4y]] (45) 

and 

> LDop(eq) : 



w -► {{Ax + 2y + 2yx + 2y 2 - 2y 3 ))(-f- w) + (-1 + x + y + 3y 2 ){4~ w) (46) 

dx dy 



14 



6 Performance 



As previously mentioned, we have used rational ODEs to exemplify the usage 
of the package commands, we hope to have clarified that usage. Now we are going 
to further exemplify the package but with a different approach. We will not be 
concerned with the input /output of the commands. In order to analyze the contri- 
bution of the present implementation, we will comment on how the package extends 
the solving capabilities of the existing solvers and how it does that. We will divide 
that into two main groups and present examples of those: 

6.1 First order ODEs with Liouvillian solutions 

The first one is the most obvious. The present implementation tackles (semi- 
algorithmically) rational first order ODEs with Liouvillian solutions and those were 
not treated neither on the original PS-approach nor in its extensions. So, this imple- 
mentation extends the applicability of these Darboux-type approaches to an wider 
range (not only on the Maple environment). 

An example of a rational first order ODE belonging to this class is: 



As previously mentioned, this equation would escape the standard PS procedure. 
Furthermore, this example (which is just an example of an infinite class) also elude 
the dsolve command. 

We would like to point out that the Maple solver can deal with some members 
of this class, for instance the Abel's ODEs but not in the general case. Our previous 
implementation extending the Prelle-Singer method [15] also solve some ODEs but 
with heuristic approaches. 

6.2 First order ODEs with elementary solutions 

The second group talks about a bonus: The present implementation solves ODEs, 
with elementary solutions, that are missed by the powerful solvers on the Maple 
commercial algebraic package and in some of the (Maple) implementations of the 
PS-approach and its extensions already mentioned. We call that a bonus since the 
present package was designed to dwell in the previously uncharted region (at least 
semi-algorithmically) of the rational first order ODEs with Liouvillian solutions but 
it fared well elsewhere also, as we shall see. But, before we go on with the example, 



dy _ y 2 (-2y + 1 - 2x + x 2 y) 
dx x 2 (— 2x + 1 — 2y + xy 2 ) 

Lsolver finds the following integrating factors for this ODE: 



(47) 




which, in turn, lead to the following solution: 




dx = C 
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let us answer the following: How can an implementation of the PS-approach [15] 
not solve ODEs with elementary solutions as we are claiming here? 

Of course, theoretically, the PS-approach would eventually catch the integrating 
factor, given enough time and computational power. Here it is where the situation 
gets complicated (at least for ODEs like the one we are about to display). In prac- 
tice, the Darboux polynomials can be of a too high a degree and that puts matters 
outside the nowadays computational capabilities (both on time and in computer 
memory). The attentive reader might ask: Why does this practical limitation not 
affect the extension here presented? The point is that the general form of the Inte- 
grating factor (11) can allow for a solution with lower Darboux polynomial degree 
to be found and that would turn the solving of the ODE feasible. Let us present 
then an example of these ODEs to clarify these ideas: 

^ = (-Ux - Uy - 28 x 3 + 14 y 3 + 40 x 4 - 58 x 5 - 19 x 2 y + 30 x 3 y+ 
ax 

-23 xV + 26 x V + 14 xy 3 + 21 x 4 y) / (x (7 x 2 + 7 x 3 + 7 x+ 

+7y + 7xy + 7y 2 + I3x 2 y + 7xy 2 + 13x 3 y + 7x 4 )) (48) 
Our program finds the following integrating factor for the above ODE: 

R = e ^ (x + 1)~ 8 (49) 
which, in turn, lead to the following solution: 

_ e^x 3 (7x + 7x 3 + x 4 - x 2 y - 7 x 2 + 7y 2 + 7 y) . . 

~ (21x 5 + x 7 + 7x 6 + 35x 4 + 35x 3 + 21x 2 + 7x + l) [ ' 

From (50) and (49), one can see that (48) admits an integrating factor R of the 
form: 

— 21x 5 + x 7 + 7x 6 + 35x 4 + 35x 3 + 21x 2 + 7x + l . , 

it = § (51) 

(x + 1) x 3 (7x + 7x 3 + x 4 — x 2 y — 7 x 2 + 7 y 2 + 7 y) 

So, in principle, the PS-approach should find it. As we pointed out above, in 
practice this does not occur since (please see (51)) one would need to find Darboux 
polynomials up to degree 7. In our case, please look at (49), we needed to use 
Darboux polynomials of degree 1. 

We remind the reader that the dsolve command contains many heuristic meth- 
ods based on Lie theory for solving ordinary differential equations [16, 17]. But, 
because of their pattern-matching based nature, there will always be an infinity of 
cases (of which the equation above is an example), for which the symmetries can 
not be found. 



6.3 Time Performance 

In order to provide information regarding the computational time performance of 
the package in solving ODEs, we present the table 1 with the ODEs we have analysed 
in the paper. The reason for this choice is that, no matter how large the sample we 
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chose to display, it will never be exhaustive, there are an infinity of possibilities. So, 
we have decided to use the ODES to which the readers have more familiarity (since 
they have seen the equations being discussed) in order for them to better analyze 
which feature influences the most the time of solution. Only one ODE analyzed in 
the paper was left out of the table, equation (48). That is so, since it is too big 
for it to fit the table in a reasonable way. The time the program took to solve this 
ODE was 1.713 s (please see the data related to the other ODEs on the table). As 
can be seen from the comparison of these data, the time the program took to solve 
this ODE is not the biggest of them all, even the ODE itself being the "scariest" . 
This shows that not only the size of the ODE matters to the time of solving it. All 
these times were obtained on a machine with an Intel Celeron M processor 360 (1.40 
GHz), 256 MB DDR SDRAM, running Windows XP Home. 



ODE 


Solution 


Time 


dy !/(l+x) 
dx x—xy—y 2 -\-x 2 


(C + Ei(l,-Z±*)) (y + x) + (y-y 2 )e X -f L =0 


0.380s 


dy _ y 2 (yax+yb+c) 
dx d 2 x 2 +2axb+b 2 


j e WW*)* {ya * + f+ c) dx +Cy-0 


4.096s 


dy _ y(i+x) 

dx x~y 2 +x 2 


iy/lerf^'^*) e^y+C-0 


0.550s 


dy -l+x+y+3y 2 
dx 2(2 x+y+xy+y 2 —y 3 ) 


(-l + 2y)' 2 / \ 

j x + y 2 e 4( x+y *) lyfRerf - x +2lL C - 

\2y/x+y 2 J 


1.002s 


dy _. y 2 (-2y+l-2x+x 2 y) 


j c ^(- 2 ^- 2 ^)^ | c-0 


1.562s 


dx x 2 {~2x+l-2y+y 2 x) 





Table 1: ODEs analyzed in the paper, their solutions and the corresponding time the 
program takes to solve them. 



7 Conclusions 

The Prelle-Singer approach is a semi-decision algorithm that can be used to solve 
analytically rational first order ordinary differential equations. Some implementa- 
tions of it have been produced [18, 15]. In [15], we have introduced an extension 
to the PS-approach and we have presented an implementation of a package that 
extends the applicability of the PS-type methods in order to tackle some ratio- 
nal first order ordinary differential equations whose solutions contain a sub-class of 
(non-elementary) Liouvillian functions (therefore, outside the scope of the original 
PS-approach). 

Here, we have implemented an extension presented on [1]. This extension now 
covers (in the same semi-decision way as the original PS approach does but for 
a wider class of differential equations) all rational first order ordinary differential 
equations with Liouvillian solutions. So, this enlarges the scope of the PS-type 
procedures previously implemented [18, 15] to allow it to tackle rational first order 
ODEs with Liouvillian solutions (see subsection 6.1). But, furthermore, it is worth 
mentioning that, as an additional bonus to this, we have enhanced the ODE solv- 
ing capabilities in Maple, even for the cases where the ODE presents elementary 
solutions (a sub set of the Liouvillian ones). Since, in the present implementation, 
the integrating factor has a different general form: it is written with an exponential 
part and, in some cases, that translates to computational speed in the determination 
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of it since the integrating factor can be found with Darboux polynomials of much 
lower degree and that can be the difference between failure and success in solving 
the rational first order ordinary differential equations (see subsection 6.2). 

In other words, apart from the theoretical extension that the present package 
now makes available for the user, thus making the analytical solving of a new class 
of rational first order ordinary differential equations possible, it also improves the 
Maple solving capabilities in the arena where the dsolve Maple command (even with 
some of its non-default enhancements turned on, e.g. the Lie package) and some 
previously made implementations of the PS-approach and its extensions [18, 15] 
already act. 

As further development for this work, we are currently working on extensions to 
the package to tackle ODEs containing elementary functions (in the same fashion 
Shtokhamer [6] used for the usual PS procedure). 
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